High-throughput transcriptome sequencing reveals the key stages of cardiovascular development in zebrafish embryos

Background The cardiovascular developmental process is a tightly regulated network involving multiple genes. The current understanding of the molecular mechanism behind cardiovascular development is insufficient and requires further research. Results Transcriptome sequencing of three developmental stages in zebrafish embryos was performed and revealed three key cardiovascular developmental stages. Then, the differentially expressed genes (DEGs) involved in cardiovascular development were screened out. The three developmental stages were 18 (T1), 24 (T2), and 42 h post fertilization (hpf) (T3), and the three stages were confirmed by detecting differences in expression between cardiomyocyte and endothelial marker genes (cmlc2, fli1) using in situ hybridization, which represents the characteristics of cardiovascular development. Thousands of DEGs were identified using transcriptome analysis. Of them, 2605 DEGs were in T1-vs-T2, including 2003 up-regulated and 602 down-regulated genes, 6446 DEGs were in T1-vs-T3, consisting of 4608 up-regulated and 1838 down-regulated genes, and 3275 DEGs were in T2-vs-T3, including 2420 up-regulated and 855 down-regulated genes. There were 644 common DEGs and 167 common five-fold higher differentially expressed genes (HDEGs) identified, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID). Significant differences was observed in the levels of gene expression among different developmental stages in multiple GO terms and KEGG pathways, such as cell migration to the midline involved in heart development, cardiovascular system development, circulatory system process for biological processes of GO terms; and cardiac muscle contraction, adrenergic signaling in cardiomyocytes for KEGG pathways. These results demonstrated that these three stages were important period for the development of the cardiovascular system. Lastly, we used quantitative real-time PCR (qPCR) to validate the reliability of RNA-sequencing by selecting 21 DEGs. Conclusions These results demonstrated that these three stages represented the important periods for cardiovascular system development of zebrafish and some candidate genes was obtained and provided a solid foundation for additional functional studies of the DEGs. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08808-x.


Background
Proper vascularization is necessary for embryonic development and adult survival. In zebrafish, the cardiovascular system develops before other organ systems during embryogenesis, with circulation beginning approximately 24 h post-fertilization (hpf) [1,2]. Cardiovascular development Open Access *Correspondence: lixiaoyu65@263.net 1 College of Life Science, Henan Normal University, Xinxiang 453007, Henan, China Full list of author information is available at the end of the article is a dynamic, spatiotemporally, and transcriptionally regulated process involving thousands of genes [3].
Zebrafish are increasingly popular in many fields of scientific research due to their rapid growth, relatively short reproductive period, and their ease of genetic manipulation [4]. They are a particularly good model animal for studying vascular development since they have transparent embryos, can be manipulated genetically [5], and their embryos can survive without blood circulation for approximately seven days post-fertilization. Therefore, it is easy to screen for vascular mutants.
In the present study, we focused on investigating important cardiovascular system-related genes during early embryogenesis, based on the zebrafish developmental model. The differentially-expressed genes (DEGs) among three key developmental stages (18,24, and 42hpf ) were analyzed using integrated transcriptome sequencing. The results demonstrated that these three stages were important for cardiovascular system development. This study provided a preliminary basis for additional exploration of DEG functions, particularly common DEGs.

Characteristics of cardiovascular development at three developmental stages
Three key developmental stages were selected, including 18, 24, and 42hpf (Fig. 1). The characteristics of cardiovascular development at three stages were demonstrated through in situ hybridization using the antisense probe of cardiomyocyte and endothelial marker genes (cmlc2 and fli1 respectively, the probe primers was listed in Supplementary Table 2). A pair of primordia were generated at 18 hpf, one on each side of the midline (Fig. 1a). The primordia on both sides fused to form a single linear heart tube at 24 hpf (Fig. 1b); this is when the heart begins to contract and blood starts to flow. By the time the embryo reaches 42 hpf, the heart tube undergoes looping morphogenesis, the heart chambers are visible, and the ventricle and atrium have been separated (Fig. 1c). The expression of endothelial marker genes (fli1) was only found in the main vessels at 18 hpf and 24 hpf embryos (Fig. 1d, e). However, when the embryo developed to 42 hpf, it was found in intersegmental vessels (ISVs), which connect the dorsal arteries and veins ( Fig. 1f) and ensure that blood flows throughout the body. These results revealed that these three stages are so important for the development of the cardiovascular system that some genes of the cardiovascular system developmental network are likely active. Therefore, comparative transcriptome analysis across three stages can provide a foundation for further research into the molecular mechanisms of genes associated with specific biological processes.

RNA sequencing and de novo transcriptome assembly
To identify gene transcription profiles of three key developmental stages, RNA sequencing was performed using Fig. 1 Expression of cardiomyocyte and endothelial marker genes (cmlc2, fli1) at three different stages of zebrafish embryo development (18, 24, 42 hpf ). a-c Expression of cmlc2 at 18hpf (a), 24hpf (b), and 42hpf (c). A pair of primordia were generated at 18 hpf, one on each side of the midline (a), The primordia on both sides fused to form a single linear heart tube at 24 hpf (b), the heart tube undergoes a loop and the heart chamber was visible (c). V, ventricle; A, atrial. d-f Expression of fli1 at 18 hpf (d), 24 hpf (e). and 42hpf (f) an Illumina HiSeq2500 system. The percentage of bases of Q20 ≥ 20 reached 99% (Fig. 2c). The correlation coefficient between the two replicates exceeded 95% (Fig. 2a), demonstrating sufficient repeatability between the two parallel experiments. The 80%-100% gene coverage accounted for approximately 70% of total genes (Fig. 2b). These results demonstrate that the sequencing data was credible enough for further analysis.

Identification of DEGs among three key developmental stages
Thousands of DEGs were found to be up-regulated or down-regulated according to the result of transcriptome analysis across three different developmental stages. (Supplementary table 3 , 4, 5). There were 2605 DEGs in T1-vs-T2, including 2003 up-regulated genes and 602 down-regulated genes, 3275 in T2-vs-T3, including 2420 up-regulated genes, and 855 down-regulated genes, and 6446 in T1-vs-T3, including 4608 up-regulated genes and 1838 down-regulated genes (Fig. 3a). T1-vs-T3 had the most DEGs compared to T1-vs-T2 and T2-vs-T3, while T1-vs-T2 had the fewest. Interestingly, there were more up-regulated DEGs than down-regulated DEGs in each comparative group. The number of DEGs up-regulated in T1-vs-T2 was 3.3 times higher than the number of down-regulated ones, 2.8 times higher in T2-vs-T3, and 2.5 times higher in T1-vs-T3. Volcano plots indicated that the number of up-regulated DEGs was significantly higher than the number of down-regulated DEGs (Fig. 3).

GO term and KEGG pathway enrichment analysis of DEGs
To further understand DEGs function during three different developmental stages, all DEGs were subjected to GO term analysis and KEGG pathway enrichment analysis [6] (Fig. 4). GO term analysis demonstrated that DEGs in T1-VS-T2 were primarily enriched in some biological processes (p < 0.05), two of them were cell migration to the midline involved in heart development (GO 0,003,318), cardioblast migration to the midline involved in heart field formation (GO 0,060,975) (Fig. 4a, Supplementary table 6). The results indicated that a pair of primordia on both sides migrated to the midline and fused to form an intermediate cone when embryos developed from 18 to 24hpf, which were consistent with the results of in situ hybridization of cardiomyocyte marker gene cmlc2 (Fig. 1). Three of GO terms which enriched in T1-VS-T3 were circulatory system process (GO 0,003,013), cardiovascular system development (GO 0,072,358), and circulation system development (GO 0,072,359) ( Fig. 4c  and formed the chamber when embryos developed from 24 to 42hpf. The results of KEGG analysis indicated that the DEGs across three groups were enriched in some pathways related to the circulatory system, such as cardiac muscle contraction, adrenergic signaling in cardiomyocytes and PI3K-AKT signaling pathway ( Fig. 4 Supplementary table 9, 10, 11). The results demonstrated that the three stages were important for the cardiovascular system development, which provide the basis for screening the genes involved in cardiovascular system development.

Identification of common DEGs in three stages
To further understand the DEG function in the three stages, common DEGs and common HDEGs of the three comparative groups were identified (Fig. 5a, b). There were 644 common DEGs (Fig. 5b); of them, 568 genes were up-regulated and 76 genes were down-regulated. There were 167 common HDEGs (Fig. 5a), including 161 up-regulated genes and 6 down-regulated genes. Interestingly, hierarchical clustering analysis of the DEGs based on normalized FPKM values indicated that common DEGs (Fig. 5c) and common HDEGs (Fig. 5d) were primarily up-regulated, while only a few were downregulated. Additionally, differences in T1 vs T2 were less important than differences in T1 vs T3 and T2 vs T3.To further explore the role of DEGs in the three comparative groups, we performed KEGG pathway enrichment analysis on 644 common DEGs and 167 common HDEGs. The results of GO analysis demonstrated that 644 common DEGs were enriched in some terms related to heart development, such as heart process, cardiovascular system development, circulatory system development, circulatory system process (Fig. 5e, Supplementary table 12), while 167 common HDEGs in cardiac myofibril assembly, cardiac cell development, cardiac muscle cell development, heart process and cardiac muscle cell differentiation (Fig. 5g, Supplementary table 13). The results of KEGG enrichment analysis demonstrated both 644 common DEGs and 167 HDEGs were enriched the cardiovascular-related pathways, such as cardiac muscle contraction, adrenergic signaling in cardiomyocytes and PI3K-AKT signaling pathways (Fig. 5f, h, Supplementary table 14,15). In conclusion, these three developmental stages were possibly crucial for cardiovascular development.

Identification of cardiovascular-related DEGs
Expression profiles demonstrated that common DEGs were associated with multiple pathways. While common DEGs and HDEGs were most significantly enriched in the circulatory system, they were also involved in other pathways, such as protein digestion and absorption and adrenergic signal transduction in cardiomyocytes. Therefore, we selected the genes enriched in the circulatory system and cardiovascular disease for hierarchical cluster analysis (Fig. S1), and the results showed that there were 42 DEGs listed in Table 1 and 11 HDEGs listed in  Table 2.

Characterization of DEGs expression in different tissues
To further investigate the potential function of certain genes during cardiovascular development, their expression in different tissues was characterized. Most of them were highly expressed in the heart or muscle. For example, the expression level of cacnb4 was highest in the brain, followed by the eyes; the expression of Slc8a1a was highest in the heart, followed by the eyes; the expression of desma was highest in the gills, followed by the heart; the expression of lama2 was highest in the muscles, followed by the heart; the expression of ryr2b was highest in the heart, followed by the muscles and eyes; and the expression of has2 was highest in the muscles, followed by the eyes and the heart (Fig. 7). The results suggest that these genes likely play a key role in cardiovascular development.

Analysis of the spatial expression of DEGs
According to the results of qPCR verification in three stages, most genes showed an increasing trend across the three stages. The results of expression in different tissues found that ryr2b, lama2, scl8a1a, has2, and desma had higher levels of expression in the muscle or the heart than in other tissues. According to the annotation of the ZFIN database, cntn2 was weakly expressed in the intersegmental vessels, so we chose these 6 genes to detect the spatial expression by the whole embryo in situ hybridization. The results demonstrated that cntn2 expressed weakly in intersegmental vessels and strongly expressed in nerve tissues (Fig. S3). Desma was expressed in myotomes, hapaxial muscles, pectoral fin muscles, head muscles, and the heart (Fig. S2). Has2 was mainly expressed in the heart, posterior somites, and tail bud (Fig. S4). Lama2 was mainly expressed in the somites (Fig. S5), Ryr2b was mainly expressed in the heart (Fig. S6), and Slc8a1a was mainly expressed in the heart (Fig. S7).

Transcriptome analysis in zebrafish embryogenesis
The rapid development of high throughput sequencing technology has allowed us to better understand the mechanism of various biological processes, including those of zebrafish. Previous transcriptome analysis primarily focused on early developmental stages [4,[7][8][9]. Results demonstrated that stage-specific genes are typically preferentially expressed during these stages and that these genes could be used as molecular markers to distinguish one stage from another [9]. This indicates that stage-preferential genes could have close relationships with biological processes. Since the morphological changes of the heart and blood vessels are significant at 18 hpf, 24 hpf, and 42 hpf in zebrafish embryos, the genes preferentially expressed during these three stages must be related to the development of heart and blood vessels. The aim of this study was to screen out genes related to cardiovascular development through transcriptome sequencing of these three stages.

Important genes involved in cardiovascular development were identified
Transcriptome analysis of three developmental stages allowed us to screen out 42 genes that could be involved in cardiovascular development. Of these genes, atp1a2a regulates zebrafish right-left patterning [10] and plays a key role in establishing proper cardiac laterality in zebrafish [11]. Additionally, certain genes screened through the transcriptome analysis are involved in the cardiovascular developmental process, including cardiac contractility [11][12][13][14], cardiomyocyte differentiation [14], vascular stability [15], hematopoiesis [16], PAAs (The pharyngeal arch arteries) development [17] and heart development [18]. Some genes are also expressed in the heart [17,19], but their functions are largely unknown and require additional study.
In this study, we used whole-mount in situ hybridization to identify the expression of 6 genes in zebrafish embryos. Of these, desmin was expressed in the head muscles, trunk muscles, pectoral fin muscles, in developing hearts and fins. These results are consistent with previous reports [20,21]. Desmin has a muscle-specific expression pattern, and is expressed in smooth muscles, skeletal muscles, and myocardium [20]. However, its function in developing hearts is still unclear and requires further research. We found that cntn2 was mainly primarily expressed in neurons and weakly expressed in intersegmental blood vessels. These results are consistent with those of previous studies [21]. Other research has reported that cntn2 is involved in branchiomotor neuron migration [22][23][24], spinal cord regeneration [25], and axon fasciculation [26]. Nevertheless, the role of cntn2 in intersegmental blood vessels has yet not been reported. Since it is expressed in intersegmental blood vessels, it could play a role in the development of intersegmental blood vessels, but this requires additional study. We determined that lama2 was primarily expressed in the trunk muscles, which was consistent with previous studies [27,28]. Laminins are the main components of basement members and are required for the development of embryonic and adult tissues [29]. It consists of three chains (alpha, beta, and gamma) which are encoded by different genes. They have different tissue-specific expression patterns and have different functions, which are described in previous reports [28]. Existing reports have demonstrated that laminins plays several roles in adults and embryos, including blood vessel formation [30], notochord formation [31], neuron migration [32], and neurogenesis [33,34]. Interestingly, laminin alpha chains have overlapping roles in the formation of notochords and blood vessels [30]. However, little is known about the role of lama2 in blood vessel formation. Therefore, this requires further study. In addition, our results also showed that ryr2b was dominantly expressed in the heart, which is consistent with previous studies [35]. However, its function in the heart has not yet been reported, and further research is necessary. Additionally, results of the heart-specific expression patterns of two genes detected in this study (slc8a1a and has2) was consistent with previous reports in which slc8a1a was reported to be required for proper cardiac morphogenesis, while Has2 plays a key role in vasculature formation [36]. In conclusion, we screened out the DEGs from three key developmental stages of zebrafish embryonic cardiovascular system using transcriptome analysis, some genes, including slc8a1a and has2, were previously reported to be involved in cardiovascular development. Some genes were expressed in cardiovascular development, but their function is unknown and need more study in the future. To this end, our research provides more candidate genes for future work and to better understand the role of these genes in embryonic cardiovascular development.

Animal husbandry and embryo collection
Zebrafish (Danio rerio) were purchased from the Institute of Hydrobiology, Chinese Academy of Science (Wuhan, China). Wild-type (AB line) zebrafish were bred and maintained according to standard procedures [37] and staged as described [38]. All animal processing was approved by the Institutional Animal Care and Use Committee of Henan Normal University.

RNA extraction, library construction, and sequencing
Total RNA was isolated using a Trizol reagent kit (Cwbiotech, Beijing, China) according to the manufacturer's instructions. The resulting RNA was further sequenced using Illumina HiSeq2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Data filtering and de novo assembly
Raw reads containing adapters or low-quality bases produced from the sequencing machines were first filtered using fastp [39]. Reads mapped to the ribosome RNA (rRNA) database using the short reads alignment tool Bowtie2 [40] (version 2.2.8) were removed. The remaining reads were mapped to the reference genome using HISAT2. 2.4 [41]. The mapped reads from each sample were assembled using StringTie v1.3.1 [42].

Quantitative real-time PCR (qRT-PCR) validation
To validate the reliability of the transcriptome sequencing, qPCR was performed for 21 DEGs. The primers were designed using Primer 5.0 and sequences were provided in Supplementary Table 1. Total RNA was extracted from 50 embryos or tissuse from two adult zebrafish. First strand cDNA was reversely transcripted from 1 μg RNA Fig. 6 Validation of DEGs by qPCR. The expression of 18 hpf was set as 1, and the relative expression level of the other two stages (24 hpf and 42 hpf ) was calculated for several genes. The asterisk indicates significantly different expression levels (*P < 0.05; **,*** and ****P < 0.01) with HIFIScript 1st strand cDNA synthesis kit (Cwbiotich, China). The independent experiment was performed in triplicate. The rpl13a gene was used as the internal control. The reaction system (10 μl) comprised 0.2 μl cDNA,5 μl SYBR green master mix, 0.4 μl each primer (20 μM), and 4 μl nuclease-free water. The PCR conditions were as follows: 10 min at 95 °C, followed by 45 cycles at 95 °C for 15 s, 60 °C for 60 s, and a cooling stage at 4 °C. The 2 −ΔΔCt method was used to analyze expression levels.

Expression profile of DEGs in different tissues
To characterize the expression profile of DEGs in different tissues, qPCR was performed and used eight tissues from adult zebrafish as templates. The reaction system and PCR conditions were the same as those above.

Spatial expression of DEGs
To further verify their spatial expression, six genes likely involved in cardiovascular development were selected to determine their expression sites in embryo development via whole-mount in situ hybridization. The primers of the probe were presented in Supplementary Table 2.
The PCR fragments were cloned into a pGEM-T easy vector and sequenced to verify whether the fragment was ligated to the vector. The antisense and sense probe were transcripted by SP6 or T7 to generate the DIG-labeled probe using the linearized vector, including PCR target fragments as a template. Whole-mount in situ hybridization was performed as previously described [44].

Statistical analysis
SPSS 20.0 was used for statistical analysis. Data was presented as mean ± SEM. Difference between groups we analyzed using the Student's t-test. P values < 0.05 were considered to be statistically significant.